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Abstract. - We present the results of molecular dynamics simulations of SPC/E water con- 
fined in a realistic model of a silica pore. The single-particle dynamics have been studied at 
ambient temperature for different hydration levels. The confinement near the hydrophilic surface 
makes the dynamic behaviour of the liquid strongly dependent on the hydration level. Upon 
decrease of the number of water molecules in the pore we observe the onset of a slow dynamics 
due to the "cage effect" . The conventional picture of a stochastic single-particle diffusion process 
thus looses its validity. 



The normal limit of supercooling and the possibility of vitrification of water is the subject 
of a longstanding scientific debate ||l|, ||. At Th ~ 236 K homogeneous nucleation, due 
to the presence of impurities, drives bulk supercooled water into its solid, crystalline phase. 
Experiments sustain the hypothesis that the amorphous phase can be connected to the normal 
liquid phase through a reversible thermodynamic path Q. Recent molecular dynamics (MD) 
simulations of bulk supercooled SPC/E Q water also showed a kinetic glass transition as 
predicted by mode coupling theory (MCT) |^ at a critical temperature Tc ^ Ts where 
Ts is the singular temperature of water which is T5 = 22%K or 49 degrees below the 
temperature of maximum density. 

A comparison of the behaviour of the bulk liquid with the same liquid in a confined environ- 
ment is highly interesting for the development of both biological and industrial applications. 
Understanding how the dynamics of liquid water is perturbed by the interaction with hy- 
drophilic or hydrophobic substrates at various levels of hydration lays the physical basis for 
the prediction of stability and enzymatic activity of biologically important macromolecules |^ . 
Many experimental studies have been performed on confined and interfacial water. They 
all showed evidence of a substantial degree of slowing down of water when confined in the 
proximity of a polar surface Q. Particularly significant in this context are the indications of 
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a transition of adsorbed water to a glassy state, which is supposedly driven by the protein 
surface . MD calculations found a reduced self diffusion coefficient of water in contact with 
solid hydrophilic surfaces, when compared to bulk water [ pi] . Moreover, a recent simulation 
study of water close to the surface of a protein evidenced a typical spectral glassy anomaly, 
the so called boson peak which the authors related to protein-solvent coupling. These 
and other studies suggest the possibility of a common underlying molecular mechanism for the 
slowing down of the single particle dynamics of interfacial water. Nevertheless, many questions 
concerning the substrate-induced perturbation of the dynamics of water are still unanswered. 
Some of these questions can be addressed by MD simulation techniques due to their power as 
a very versatile microscopic tool. 

In this letter, we present the results of MD simulations of SPC/E water confined in a 
cylindrical silica pore. The confining system has been carefully designed to mimic the pores 
of real Vycor glass. We chose this particular system, because Vycor is characterized by a 
quite sharp distribution of pore sizes with a small average diameter 40 ± 5 A), because 
the pore size does not depend on the hydration level (no swelling occurs), and because the 
surface is strongly hydrophilic. For the same reasons, this system has been studied extensively 
by a variety of experiments Here, we focus on the role of hydration. In 

particular, we try to understand the extent to which the changes in the system dynamics upon 
lowering the hydration level in a hydrophilic silica pore are equivalent to supercooling the bulk 
phase To our knowledge, this is the first time that the dynamics of water in a realistic 

pore model is calculated in such detail. The results presented here can help to shed light on 
the microscopic mechanism governing the dynamics of water close to a hydrophilic surface. 

The MD simulations started from a silica cube containing a cylindrical cavity. In the 
system preparation step a vitreous silica cubic cell was first created by melting a /3-cristobalite 
crystalline structure and then quenching the liquid down to ambient temperature [ p^ p^ . 
Inside the glass, which is formed by this process, a cylindrical cavity with a radius of i? = 20 A 
is created by appropriate removal of atoms. Next, all those remaining silicon atoms which are 
bonded to less than four oxygens are removed. With this procedure two kinds of oxygen 
atoms are obtained, nonhridging oxygens (NBOs), bonded to only one silicon ion, and bridging 
oxygens (BOs), bonded to at least two silicons. The NBOs are saturated with hydrogen 
atoms, similar to the experimental process before hydration fl^ , which leads to free hydroxyl 
groups on the pore surface. With this process, we create a cavity with a rough surface, 
whose geometry, size, and microscopic structure are typical for a pore in Vycor glass. The 
block of silica provides a rigid matrix for the simulation of SPC/E water molecules which are 
then introduced into the pore. Water interacts with the atoms of the rigid matrix by means 
of an empirical potential model ||l^, where different fractional charges are assigned to BOs 
(-0.629|e|), NBOs (-0.533|e|), silicon (+1.283|e|) and surface hydrogen atoms (-|-0.206|e|). 
Additional Lennard- Jones potential functions between the oxygen sites of water and BOs and 
NBOs complement the water- Vycor interaction [|l8| . Since the simulations are time-consuming, 
we have used the shifted force method with a cutoff at 9 A. We did check, however, that the use 
of a larger cutoff or Ewald summations docs not change the trend of the obtained results . 

There is experimental evidence that the average density of water confined in real Vycor 
glass at full hydration is 11% less than the density of bulk water at ambient conditions |2^. In 
order to match this value 2700 water molecules are needed in our pore, this number defines our 
full hydration. We have conducted simulations in the NVE ensemble at ambient temperature 
(T = 298X) for four different hydrations of the pore. Simulations with 500, 1000, 1500, 
2000 and 2600 water molecules correspond, based on the above consideration, to hydration 
levels of 19%, 37%, 56%, 74% and 96%, respectively. The timestep for the integration of the 
molecular trajectories is 1 fs. Thermal equilibrium was achieved through the use of a Berendsen 
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thermostat 1 21 1 ; the equihbration process has been carefully monitored via the time dependence 
of the potential energy. Further details of the simulation can be found in ref . , where static 
properties of the system are studied, and where it is shown that the interaction model used here 
leads to satisfactory agreement with the experimental jlj] site-site pair correlation functions. 

We now discuss some of the most significant single particle dynamical properties of water 
from our simulations for different hydration levels at ambient temperature. In Fig. 1 we show 
the mean square displacement (MSD) of the center of mass motion along the non-confined 
z-direction. After an initial ballistic diffusion, water molecules in the bulk phase at ambient 
temperature (not shown) enter the Brownian diffusive regime (characterized by a slope of 1 
in the log- log plot). Our data show at intermediate times an obvious flattening of the MSD 
curve, which increases with decreasing hydration level. This is the signature of the so-called 
cage effect in the system, which describes the trapping of molecules in the cage formed by 
their nearest neighbors at intermediate times after the initial ballistic regime. Once the cage 
has relaxed, molecules enter the normal diffusive regime. In the inset of Fig. 1 we report the 
self diffusion coefficient, D, as calculated from the slope of the MSD in the diffusive regime. 
A substantial decrease of D with decreasing hydration level is evident. All values of D are 
substantially lower than the one of SPC/E bulk water at the same temperature. 

Figure 2 shows the single particle intermediate scattering function (ISF), Fs{Q,t), calcu- 
lated for the center of mass motion at the maximum of the oxygen-oxygen structure factor for 
all hydration levels investigated. As the hydration is lowered, the ISF displays an increasingly 
pronounced shoulder around 1 ps, which is the characteristic of a two-step relaxation process. 
The long time tail is highly non-exponential for all hydration levels We are able to fit 

a combination of the Kohlrausch- William- Watts (KWW) relaxation function plus a gaussian 
term 

FsiQ, t)^[l- AiQ)] exp [- (i/r,)'] + A{Q) exp [- {t/nf] (1) 

only to the ISF at the lowest hydration level. A{Q) — e"" ' ^'^ is the Debye- Waller factor 
arising from the cage effect, a is the amplitude of the vibration of the molecule inside the 
cage created by the potential barrier of the nearest neighbours. Ts and ti are, respectively, the 
short and the long relaxation times. The fitted curve is shown as the dashed line in Fig. 2. 
For bulk water, this function fitted the center of mass ISF for several temperatures For 
the shown fit we obtain a cage radius a ~ 0.44 A, which is similar to the radius obtained 
for bulk supercooled water, a ~ 0.5 A. The lower radius obtained for confined water may 
be due to the slightly higher density of water close to the surface . The short relaxation 
time Ts ~ 0.14 ps is again comparable to the bulk value ~ 0.2 ps. The long time tail is 
characterized by /3 = 0.35 and r; = 356 ps. 

For the higher hydration levels we have calculated the ISF separately for the first two layers 
of molecules close to the surface and for the remaining molecules. This choice is based on 
the shape of the density profile, which shows a double layer of water molecules close to the 
surface with density higher than the average |]l9| . This double layer extends up to 5 A from 
the surface. The inset of Fig. 2 shows, for Nw — 2600 (96% hydration), the two contributions. 
While the molecules in the adsorbed layer relax very slowly, the decay of the ISF of the 
remaining molecules is much faster and very well described by eqJ3 (bottom curves). From the 
fit we extract f3 = 0.71, r; ~ 0.5 ps, Tg — 0.17 ps, and a = 0.54 A. While the a, Ts and t; are 
similar to the values found for bulk water at ambient conditions, the (5 differs from the value, 
1, of the bulk. It turns out that the behaviour of molecules belonging to the first hydration 
layers change as a complete surface coverage is achieved. In fact for Nw = 500 the surface 
coverage is not complete and patches of water molecules are visible along the pore surface. A 
stretched exponential function is able to account also for the dynamics of clusters in a frozen 



4 



EUROPHYSICS LETTERS 



environment 1). In Fig. 3 FsiQ, t) at 19 % hydration is shown for several values of Q ranging 
from 0.5 to 4 A~ . This fi gure very clearly demonstrates the overshooting of the ISF around 
1 ps, which was also visible in Fig. 2. 

The maximum can be observed both in the z direction along the pore axis and in the xy 
plane, as is visible in the inset. The intermediate maximum of the ISF has been related to 
the so-called boson peak (BP) The BP is an excess of vibrational modes present in many 
glasses at frequencies around 1 THz. When this glassy anomaly appears in a liquid phase, it 
is usually considered as a precursor to the actual glass transition. We also performed a shell 
analysis of the dynamic behaviour at the higher hydration levels, and we found in these cases 
that the contribution to the BP in our simulation comes only from water molecules which 
are not in the first layer. This fact, and the fact that our substrate is a rigid framework, is 
an indication that the BP is a feature of liquid water, which is not induced by the substrate 
dynamics. 

Experimental signatures of our observations were detected in some confined hydrogen 
bonded complex liquids p^ . An analysis of quasi-elastic neutron scattering data of the wa- 
ter / Vycor system is consistent with a highly non-exponential relaxation behaviour . The 
BP for the water / Vycor system has been detected recently at energies around 3.5 meV [ p5[ . 

In summary, we presented MD results concerning the single particle dynamics of liquid water 
confined in a silica pore. We found evidence of glassy behaviour already at room temperature. 
On lowering the hydration level of water inside the pore, the MSD flattens at intermediate 
times due to a cage effect. Correspondingly, the ISF displays a two step relaxation behaviour 
with a highly non-exponential slow relaxation. Such behaviour is typical for a glass forming 
liquid approaching its glass transition point. At the lowest hydration level, a KWW function 
could be fitted to the ISF. In this particular system, all phenomena are strongly influenced 
by the substrate. Specifically the interaction with the hydrophilic surface seems to drive the 
liquid closer to the glass transition point. This interaction leads to a significant slowing down 
of water molecules close to the substrate. The strong distortion of the hydrogen bond network 
close to the surface will lead to a strong variation of the cage structure from one molecule to 
the other, which, in contrast to bulk water, might mask the simple behaviour suggested by 
mode coupling theory. The dynamic results moreover support the notion of two quite distinct 
subsets of molecules. One is bound directly to the substrate surface, the other consists of 
the remaining water molecules in the pore center. Near the substrate surface, the local water 
density is slightly higher |p^ , and the dynamics is severely slowed down at all hydration levels, 
whereas the retardation of the slow relaxation process is much less pronounced for the inner 
water shells, whose contribution to the slow relaxation in the ISF could be fitted to a KWW 
law. The shouldering of the relaxation law is indeed determined by the propagation of the 
perturbation of the substrate to the water molecules inside the pore. With the increase of the 
hydration level the effect is overthrown by the increasing number of molecules located far from 
the surface. 

The confinement also leads to the formation of an overshooting related to the BP at the 
lowest hydration level. Since the BP is found typically in strong glass formers |2^, the 
appearance of the intermediate maximum in the ISF at the lowest hydration would favor 
the hypothesis of a conversion to a more strong glass-forming behaviour of water as the level 
of hydration is lowered. 

P.G. wishes to thank S.-H. Chen for all the interesting discussions on this subject. 
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Fig. 1. - Log- log representation of the mean square displacement as a function of time for the 
hydration levels studied. From bottom to top the curves correspond to increasing hydration 
levels. The inset shows self diffusion coefficients, D, vs. number of molecules in the pore (filled 
circles). The line is a guide for the eye, and the filled diamond on the y-axis is the value of D 
for SPC/E bulk water at ambient conditions. 

Fig. 2. - Intermediate scattering function (ISP) of the center of mass motion at Q = 2.25 A ^ , 
corresponding to the peak of the structure factor. The ISF is shown over four decades; full 
lines correspond, from top to bottom, to increasing hydration levels. The dashed line is the 
fit of Eq. |l| to the data at 19 % hydration. From the fit we extract (3 — 0.35, ti — 356 ps, 
Ts = 0.14 ps, and a = 0.44 A (see also text). The inset shows the ISF at 96 % hydration. The 
top curve (labeled S) is the contribution to the total ISF from the shell of molecules closest 
to the substrate (15 < r < 20 A), the bottom curve (labeled C) the contribution from the 
molecules in the remaining shells (0 < r < 15 A), and the central one (labeled T) is the total 
ISF. The radius of the pore is i? = 20 A and the zero is taken at the center of the pore. The 
dashed line is the fit of Eq. |l|. From the fit we extract f3 — 0.71, t; = 0.5 ps, — 0.17 ps, and 
a = 0.54 A. 

Fig. 3. - Self contribution to the intermediate scattering function (ISF) at 19 % hydration 
for various values of the scattering vector Q in the z direction. From top to bottom, Q = 

0. 5.1.1.2.1.5.1.8.2.2.25.2.5.3.3.5.4 A-i. The inset shows the ISF for the components of Q 
along the z direction (full line), in the xy plane (dotted-dashed line), and the sum of the two 
(dashed line) for two Q values. The first three curves on the top correspond to Q — 2.25 A~^, 

1. e. the peak of the structure factor, and the three bottom ones to Q = 4 A~^. The maximum 
related to the boson peak (BP) is evident at around 1 ps. 
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